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In  this  work,  a  multi  scale  modeling  approach  has  been  developed  to  simulate  the  impact  of  woven 
fabrics  using  a  finite  element  (FE)  analysis.  A  yarn  level  of  resolution  is  used  in  the  model.  This  approach, 
referred  to  as  the  hybrid  element  analysis  (HEA)  is  based  on  decreasing  the  complexity  of  the  finite 
element  model  with  distance  away  from  the  impact  zone  based  on  the  multiscale  nature  of  the  fabric 
architecture  and  the  physics  of  the  impact  event.  Solid  elements  are  used  to  discretize  the  yarns  around 
the  impact  region,  which  transition  to  shell  elements  in  the  surrounding  region.  A  new  method  for 
modeling  the  shell  yarns  is  incorporated  that  more  accurately  represents  the  contours  of  the  yarn  cross 
section.  Impedances  have  been  matched  across  the  solid— shell  interface  to  prevent  interfacial  reflections 
of  the  longitudinal  strain  wave.  The  HEA  method  is  validated  by  first  applying  it  to  the  FE  model  of 
a  single  yarn  for  which  an  analytical  solution  is  known.  The  HEA  method  is  then  applied  to  a  woven  fabric 
model  and  validated  by  comparing  it  against  a  baseline  model  consisting  of  yarns  discretized  using  only 
solid  elements. 

©  2010  Elsevier  Ltd.  All  rights  reserved. 


1.  Introduction 

Various  applications  have  utilized  the  effective  penetration 
resistance  and  impact  energy  dissipation  offered  by  woven  fabrics 
comprised  of  high  performance  materials  such  as  Kevlar,  Zylon,  and 
Twaron.  These  applications  include  protective  clothing  for 
personnel  extremity  protection  and  materials  for  turbine  fragment 
containment.  These  high  performance  materials  have  a  high 
strength  to  weight  ratio,  high  stiffness,  and  high  strength,  making 
them  ideal  for  these  applications.  Considerable  research  has  been 
conducted  to  understand  the  intricate  architecture  of  these  fabrics, 
as  well  as  for  capturing  their  non-linear  behavior  and  complex 
interactions  during  in-plane  loading  and  transverse  impact.  Various 
analytical,  numerical  and  experimental  approaches  have  been 
employed.  Some  of  these  approaches  and  the  factors  that  affect  the 
impact  response  of  flexible  woven  composites  have  been  outlined 
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in  the  reviews  by  Tabiei  and  Nilakantan  [1  ],  and  by  Cheeseman  and 
Bogetti  [2]. 

There  are  several  techniques  used  to  represent  fabric  behavior. 
Single  scale  modeling  techniques  include  representing  the  entire 
fabric  as  homogenized  membranes  [3-5],  explicitly  capturing  yarn 
level  architecture  [6-8],  and  capturing  filament  level  architecture 
[9,10].  Neglecting  yarn  architecture,  as  is  the  case  when  repre¬ 
senting  the  fabric  by  a  homogenous  membrane  leads  to  excellent 
computational  efficiency,  but  these  approaches  often  use  over¬ 
simplifying  assumptions.  For  example  individual  yarn  failure, 
capturing  phenomena  such  as  yarn  pullout  and  sliding,  accounting 
for  transverse  yarn  compression  and  thickness  changes,  and 
computing  the  energy  dissipated  by  frictional  sliding  between 
yarns  are  not  possible  with  these  models.  Frictional  effects  have 
been  shown  to  have  an  important  effect  on  the  fabric  response  to 
impact  and  researchers  have  been  investigating  methods  to  alter 
inter  yarn  frictional  properties  with  a  view  to  improving  the  energy 
dissipation  capabilities  of  these  fabrics  [11,12].  On  the  other  hand, 
models  that  explicitly  capture  yarn  architecture  are  more  accurate 
and  can  account  for  individual  yarn  motion  and  interactions. 
However  these  models  usually  require  enormous  computational 
resources  in  order  to  be  able  to  model  multi  layer  fabric  systems  of 
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realistic  dimensions.  Such  computational  resources  are  currently 
unavailable.  Some  finite  element  (FE)  codes  also  tend  to  become 
unstable  when  dealing  with  the  enormous  number  of  degrees  of 
freedom  associated  with  such  models. 

More  recently,  multi  scale  modeling  techniques  have  emerged 
that  attempt  to  balance  the  inequity  between  accuracy  and  fineness 
of  modeling  resolution  with  computational  expense.  Barauskas  and 
Abraitiene  [13,14]  used  a  mezzo-  and  macro-mechanical  model  that 
used  thin  shell  elements  to  model  yarns  at  the  impact  point  and 
roughly  meshed  uniform  orthotropic  thin  shells  for  zones  far  away 
from  the  impact  point.  A  tied  interface  was  used  to  couple  both 
regions.  However  their  model  could  not  adequately  capture 
changes  in  the  yarn  thickness  and  cross  sectional  shape,  and  their 
model  required  a  manipulation  of  the  contact  algorithm’s  penalty 
coefficient  in  an  attempt  to  mimic  the  through  thickness  defor¬ 
mation  of  the  yarns.  The  selection  of  the  material  model  and  cor¬ 
responding  material  parameters,  especially  in  the  zones  remote 
from  the  impact  site  is  in  an  arbitrary  manner  and  this  is  reflected 
by  large  errors  in  the  wave  propagation  speeds  which  directly  affect 
the  accuracy  of  results,  especially  in  high  rate  deformation  condi¬ 
tions.  No  comparison  has  been  made  against  a  baseline  numerical 
model  in  terms  of  either  simulation  results  such  as  projectile 
velocity  history  and  energy  interaction  histories  or  simulation  run 
times,  although  the  purpose  of  such  a  model  is  to  reflect  savings  in 
computational  expense  while  maintaining  accuracy.  It  is  essential 
that  a  FE  model  be  predictive  in  nature  and  not  dependent  on 
arbitrary  adjustment  of  input  parameters  to  fit  experimental  or 
baseline  data.  Rao  et  al.  [15]  present  a  multi  scale  model  that  uses 
a  ‘ center-cross '  modeling  approach.  Yarns  are  discretely  modeled 
using  solid  elements  within  a  central  cross  or  local  region  which 
extends  up  to  the  fabric  boundaries,  while  a  homogenized  layer 
modeled  with  two  layers  of  solid  elements  is  used  in  the  global 
region,  which  is  the  region  outside  the  central  cross  shape.  This 
allows  for  some  degree  of  yarn  sliding  and  yarn  pullout,  especially 
in  cases  where  the  fabric  is  gripped  on  two  ends.  However  main¬ 
taining  local  yarn  architecture  modeled  using  solid  elements  all  the 
way  to  the  fabric  boundaries  may  be  impractical  in  terms  of 
computational  expense.  Their  model  uses  multiple  layers  of  solid 
elements  in  the  global  region  leading  to  a  far  more  computationally 
expensive  model  than  if  shell  elements  were  used  in  the  global 
region.  Their  multi  scale  model  compares  well  to  their  baseline 
numerical  model,  which  is  comprised  wholly  of  yarns  modeled 
using  solid  elements. 

As  such  the  field  of  multi  scale  numerical  modeling  of  the 
impact  of  flexible  woven  fabrics  is  fairly  recent,  and  much  research 
still  needs  to  be  conducted  in  order  to  improve  accuracy  and 
establish  a  consistent  modeling  methodology,  while  reducing 
computational  expense.  To  this  end,  a  multi  scale  numerical  model 
that  uses  a  new  technique  entitled  ‘ hybrid  element  analysis'  is 
developed.  The  HEA  method  is  compared  to  a  baseline  numerical 
model  where  a  fabric  has  been  completely  modeled  with  a  yarn 
level  resolution  using  solid  elements.  The  purpose  of  the  HEA 
model  is  to  accurately  reproduce  the  results  of  the  baseline 
numerical  model,  but  at  a  much  lower  computational  expense.  This 
study  introduces  the  HEA  approach  starting  with  a  single  yarn 
model  and  extending  it  to  a  fabric  model  with  a  yarn  level  archi¬ 
tecture  used  throughout.  This  is  referred  to  as  the  ‘ single  scale  HEA'. 

2.  Modeling  of  yarns  using  shell  elements 

Solid  or  hexahedral  finite  elements  are  a  common  choice  when 
modeling  yarns.  Fig.  1  displays  the  FE  mesh  of  an  uncrimped  woven 
yarn  from  a  plain  weave  fabric  modeled  using  solid  finite  elements. 
A  mesh  sensitivity  study  was  conducted  to  decide  a  suitable  mesh 
scheme  that  balances  computational  expense  with  accuracy  of 


Fig.  1.  FE  mesh  of  a  yarn  modeled  using  solid  finite  elements. 


results.  Six  elements  are  used  across  the  yarn  width  to  sufficiently 
capture  the  curved  shape  of  the  yarn  cross  section,  and  two 
elements  are  used  through  the  thickness.  Since  yarns  are  essentially 
very  thin  structures,  in  order  to  maintain  good  element  aspect 
ratios  a  large  number  of  elements  need  to  be  used  along  the  yarn 
length  so  that  the  element  dimension  along  the  length  is  the  same 
order  of  magnitude  as  the  element  dimension  across  the  width. 
This  is  also  important  in  order  to  maintain  smooth  contacting 
surfaces  between  two  orthogonal  yarns  that  slide  past  one  another. 
This  makes  using  only  solid  elements  a  computationally  intensive 
choice  when  modeling  yarns  of  very  large  lineal  dimensions. 

The  principal  idea  behind  using  a  shell  element  is  to  collapse 
a  three  dimensional  object  into  two  dimensions,  which  also  enables 
the  modeling  of  very  thin  structures.  Shell  elements  are  computa¬ 
tionally  less  intensive  than  solid  elements,  thereby  making  them 
useful  for  modeling  yarns.  Shell  elements  have  a  minimum  of  two 
integration  points  through  the  thickness  allowing  them  to  account 
for  both  bending  stiffness  and  thickness  changes  unlike  membrane 
elements  that  have  only  one  through  thickness  integration  point. 
The  thickness  change  in  a  shell  element  is  determined  from  the 
computation  of  the  through  thickness  strain,  which  in  turn  is 
computed  from  the  in-plane  normal  stresses  and  Poisson  ratios. 
Thus  the  transverse  compression  of  a  shell  element  is  not  updated 
separately  since  the  plane  stress  condition  mandates  that  the 
through  thickness  stress  is  zero,  rather  it  depends  on  the 
membrane  straining.  However  the  observed  thickness  changes  in 
yarns  are  mainly  due  to  transverse  compressions  as  opposed  to 
longitudinal  straining.  Thus  the  constant  thickness  of  a  shell 
element  is  a  potential  limitation  when  trying  to  use  shell  elements 
in  regions  where  capturing  transverse  yarn  compressions  is 
important,  such  as  directly  under  the  projectile  at  the  impact 
location.  However  shell  elements  are  easily  able  to  capture  the  axial 
tensile  deformations  of  a  yarn  modeled  with  solid  elements. 

Previous  approaches  [13,16]  to  modeling  yarns  using  shell 
elements  are  shown  in  Fig.  2a.  As  can  be  observed,  the  cross 
sectional  shape  is  not  realistically  modeled.  Each  shell  element  and 
its  corresponding  nodes  have  been  assigned  a  uniform  thickness  t-i 
to  t-iii.  There  are  thickness  jumps  at  the  shell  element  boundaries 
[16].  This  is  not  desirable  in  the  impact  of  textile  fabric  structures 
where  it  is  essential  to  accurately  capture  inter  yarn  sliding  inter¬ 
actions  which  require  realistic  representations  of  the  yarn  cross 
section.  Another  drawback  is  the  loss  of  yarn  material  at  the  regions 
of  the  thickness  jumps  leading  to  a  loss  in  yarn  mass.  Two  possi¬ 
bilities  of  conserving  the  yarn  mass  are  to  either  artificially  increase 
the  yarn  material  density  or  to  increase  the  thickness  of  the  shell 
elements.  The  former  leads  to  a  different  longitudinal  strain  wave 
velocity  and  prediction  of  longitudinal  stresses  within  the  yarn, 
while  the  latter  leads  to  inter  yarn  penetrations  at  the  yarn  cross¬ 
over  locations  in  a  fabric,  both  of  which  are  undesirable.  A  new 
modeling  approach  is  illustrated  in  Fig.  2b,  where  a  yarn  has  been 
modeled  using  the  same  number  of  shell  elements  across  the 
width,  with  non-uniform  nodal  thicknesses  tl— 14.  These  non- 
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Fig.  2.  Yarn  cross  section  using  shell  elements  with  (a)  uniform  nodal  thicknesses,  (b)  non-uniform  nodal  thicknesses,  (c)  FE  mesh  of  a  yarn  modeled  using  shell  elements. 


uniform  nodal  thicknesses  have  been  specified  for  each  shell 
element,  so  that  it  closely  approximates  the  cross  section  of  a  yarn. 
Such  a  modeling  approach  leads  to  a  more  accurate  representation 
of  the  yarn  geometry  and  mass.  The  preprocessor  DYNAYarn  [17] 
was  used  to  create  the  FE  mesh  and  automatically  assign  the  non- 
uniform  nodal  thicknesses  to  each  shell  element.  Fig.  2c  displays 
the  FE  mesh  of  a  yarn  modeled  using  shell  elements  which  is 
equivalent  to  the  yarn  modeled  using  solid  elements  in  Fig.  1. 

To  compare  this  new  method  of  modeling  yarns  with  shell 
elements  (see  Fig.  2b)  with  the  previous  method  (see  Fig.  2a),  a  high 
rate  yarn  tensile  test  is  simulated  using  the  dynamic  finite  element 
code  LS-DYNA.  Consider  the  case  of  a  50  mm  long  uncrimped  yarn 
gripped  at  one  end  and  pulled  at  the  other  end  in  the  x-direction 
with  a  constant  velocity  of  20  m/s.  The  simulation  is  run  for  50  ps 
which  results  in  a  net  displacement  of  1  mm  at  the  pulled  end.  A 
linearly  elastic  orthotropic  material  model  (Mat  #2  in  LS-DYNA)  is 
assigned  to  the  yarn  with  the  following  properties:  longitudinal 
elastic  modulus  (Ex)  of  89  GPa,  density  (p)  of  1440  kg/m3,  and  strain 
to  failure  (£f)  of  3.1%.  According  to  Ref.  [18],  an  orthotropic  elastic 
continuum  can  be  used  to  model  a  continuous  filament  yarn 
provided  the  transverse  moduli  ( Ey ,  Ez )  are  very  small  compared  to 
the  longitudinal  elastic  modulus  ( Ex ),  and  with  zero  Poisson  ratios 
and  shear  moduli  ( Gxy ,  Gyz,  Gzx).  This  approach  essentially  implies 
a  one  dimensional  response  similar  to  a  cable  or  beam  wherein  the 
longitudinal  response  is  the  most  predominant.  However  in  a  finite 
element  analysis,  assigning  the  shear  moduli  to  zero  could  lead  to 
element  hourglassing  or  zero  strain  energy  modes  and  so  a  small 
value  usually  needs  to  be  assigned  to  the  material.  For  lack  of  better 
data,  shear  moduli  (Gxy,  Gyz,  Gzx )  of  3280  MPa,  and  zero  Poisson 
ratios  have  been  selected  [19].  Here  the  x-axis  is  along  the  length 
direction,  the  y-axis  is  along  the  thickness  direction,  while  the  z- 
axis  is  along  the  width  direction.  Four  test  cases  are  studied  where 
the  yarn  is  modeled  with  (#1 )  solid  elements  which  serve  as  the 
baseline,  (#2)  our  new  method  with  shell  elements  which  is  also 
part  of  the  HEA  approach,  (#3)  previous  method  with  shell 
elements,  see  Fig.  2a,  and  (#4)  previous  method  with  shell 
elements  where  the  density  has  been  increased  to  account  for  the 
lost  yarn  mass  due  to  the  thickness  jumps.  Fig.  3  compares  the  total 
yarn  mass  in  all  four  cases.  The  yarn  masses  have  been  normalized 
with  respect  to  case  #1.  As  can  be  seen,  there  is  about  a  12%  loss  in 
mass  due  to  the  thickness  jumps  in  case  #3  which  is  rectified  in 


case  #4  by  increasing  the  yarn  density  from  1440.00  kg/m3  to 
1636.58  kg/m3.  Table  1  displays  the  yarn  thickness  at  various 
locations  across  the  yarn  width  for  each  case,  with  reference  to  Figs. 
1, 2a  and  b  respectively.  Before  we  compare  the  results  between  the 
four  cases,  we  will  first  review  some  1  -d  wave  propagation  theory 
to  formulate  some  analytical  predictions  for  comparison  with  the 
numerical  predictions,  with  respect  to  the  case  of  a  yarn  held  at  one 
end  and  pulled  at  a  constant  velocity  at  the  other  end. 

Wave  propagation  in  bars  is  well  approximated  by  the  one 
dimensional  wave  equation  [20] 


8 2_u  _  292u 
W  ~  C  Ox1 


(1) 


which  can  be  factored  into  a  first  order  equation  for  this  case 
study  as 


0 u  _  du 
dt  =  °dx 


(2) 


Yarn  Model 

Fig.  3.  Comparison  of  yarn  mass. 
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Table  1 

Thicknesses  across  the  yam  width. 


Baseline  (solids) 

Previous  method 
(shells) 

New  method  (shells) 

Node 

t  (mm) 

Element 

t  (mm) 

Node 

t  (mm) 

tl 

0.049366 

t-i 

0.049366 

tl 

0.049366 

t2 

0.084034 

t-ii 

0.084034 

t2 

0.084034 

t3 

0.106978 

t-iii 

0.106978 

t3 

0.106978 

t4 

0.115000 

t4 

0.115000 

where  V  is  the  displacement  in  the  x-direction,  and  ‘c’  is  longitu¬ 
dinal  strain  wave  speed  given  by 


c 


(3) 


Here  ‘Fx’  is  the  longitudinal  tensile  modulus  and  y  is  the 
material  density.  The  longitudinal  strain  in  the  yarn  is  given  by 


(4) 


where  the  prescribed  displacement  ‘u’  at  the  pulled  end  can  be 
described  as  a  function  of  time  by  a  generalized  polynomial  as 


u(t)  =  at +  bt2  +  ct3  +  ...  (5) 

Then  for  the  case  of  a  prescribed  linear  displacement  function 
with  respect  to  time  (constant  velocity),  we  can  obtain  the  longi¬ 
tudinal  strain  from  Eqs.  (2),  (4)  and  (5)  as 


(6) 


where  ‘a’  is  the  constant  velocity.  The  reaction  force  at  the  clamped 
yarn  end  is  given  by 


2 EXA  du 

Fx  =  —  di  (7) 

where  ‘A’  is  the  cross  sectional  area  of  the  yarn.  The  reaction  force 
can  also  be  calculated  from  the  impulse-momentum  equation 
which  in  its  generalized  form  is  given  by 

=  FAt  (8) 

where  force  ‘F  acting  over  a  short  time  duration  ‘At’  causes 
a  reduction  in  the  velocity  ‘v’  of  travelling  mass  ‘m’  from  the 
previous  time  step  ‘n  -  1’  to  the  current  step  ‘n’.  Through  a  simple 
derivation,  Eq.  (8)  leads  to  an  alternate  form  for  the  reaction  force  at 
the  clamped  yarn  end  of  Eq.  (7)  as 


Fx  =  2  pAc2ex 


(9) 


The  factor  of  ‘2’  used  in  Eqs.  (7)  and  (9)  accounts  for  the  arrival 
and  almost  immediate  reflection  of  the  longitudinal  strain  wave  at 
the  clamped  end.  If  the  yarn  is  considered  to  be  composed  of  an 
infinite  chain  of  elements,  each  element  will  show  a  sharp  jump  in 
their  stress  level  every  time  the  longitudinal  strain  wave  passes  it. 
The  extent  of  the  stress  increment  depends  on  the  nature  of  the 
prescribed  displacement  at  the  pulled  end.  For  the  constant 
velocity  prescribed  displacement  condition  considered  here,  the 
stress  increases  by  a  constant  increment  with  each  passing  of  the 
longitudinal  wave.  This  leads  to  a  stepped  formation  in  the  stress 
plot  which  is  clearly  observable  from  the  results.  The  calculation  of 
longitudinal  stress  at  a  particular  location  ‘x’  along  the  yarn  length 
follows  from  the  calculation  of  strain  as 


_  Ex  du 

&X  —  —  TTT 

C  dt 


(10) 


The  distance  between  the  steps  will  depend  on  the  location  of 
the  element  along  the  yarn  length  and  the  corresponding  time  it 
takes  for  one  reflection  of  the  longitudinal  wave  from  the  point  of 
interest  to  the  pulled  end,  and  also  to  the  fixed  end. 

The  internal  energy  or  elastic  strain  energy  of  the  yarn  is  given  by 

i 

SEyam  =  2J  ^x^x^dx  (H) 

0 

The  kinetic  energy  of  this  yarn,  held  at  one  end  and  pulled  at  the 
other  end  with  constant  velocity,  is  given  by 

KEyarn  =  \mv2  =  l(PAl)(cex)2  (12) 

where  T  is  the  instantaneous  yarn  length,  i.e.  the  length  of  yarn 
behind  the  front  of  the  longitudinal  strain  wave.  It  is  important  to 
note  in  the  above  equations  that  the  quantity  ex,  and  consequently 
all  other  terms  that  depend  on  it,  varies  with  time  depending  on 
the  propagation  of  the  longitudinal  wave.  Thus  Eq.  (6)  actually 
represents  the  increment  by  which  the  longitudinal  strain  keeps 
increasing  in  a  step  wise  manner  for  each  passing  of  the  longi¬ 
tudinal  wave  at  any  fixed  location  along  the  yarn  length.  Fig.  4 
compares  the  results  between  all  cases.  As  can  be  observed  the 
new  method  of  modeling  yarns  with  shell  elements  which  is  used 
within  the  HEA  approach  closely  agrees  with  the  baseline  case. 
There  is  a  good  agreement  between  the  yarn  internal  energy  for 
cases  #1  and  #2,  as  seen  from  Fig.  4a.  The  peak  yarn  kinetic  energy 
from  Fig.  4b  for  cases  #1  and  #2  closely  agrees  with  the  analytical 
prediction  obtained  using  Eq.  (12),  with  less  than  a  two  percent 
difference  when  normalized  with  respect  to  the  analytical 
prediction.  The  fluctuation  of  the  kinetic  energy  between  a  near 
zero  value  and  the  peak  is  in  accordance  with  the  propagation  of 
the  longitudinal  wave  through  the  yarn  that  reflects  at  the 
boundaries.  Using  Eq.  (7),  the  first  increment  of  the  reaction  force 
at  the  clamped  end,  viz.  at  the  arrival  and  immediate  reflection  of 
the  longitudinal  wave  is  22.1  N.  The  reaction  force  grows  by  this 
increment  with  each  reflection  of  the  longitudinal  wave  at  the 
boundary.  Comparing  the  reaction  force  predictions  of  cases  #1 
and  #2  from  Fig.  4c,  we  see  that  both  cases  closely  agree  with  the 
analytical  prediction,  with  less  than  half  a  percent  difference  when 
normalized  with  respect  to  the  analytical  prediction.  However  for 
all  the  results,  we  see  that  cases  #3  and  #4  differ  from  the  baseline 
results.  Fig.  4d  compares  the  error  in  predictions  of  cases  #3  and 
#4  with  respect  to  the  analytical  predictions.  Both  the  magnitude 
and  percentage  of  error  in  the  reaction  force  at  the  fixed  end  and 
the  time  instant  at  which  the  longitudinal  wave  reflects  from  the 
fixed  end  have  been  shown.  Because  of  the  smaller  longitudinal 
wave  speed  in  case  #4  caused  by  the  increased  yarn  density,  the 
peaks  are  shifted  to  the  right  and  lag  behind  the  other  cases. 
Consequently  the  error  in  both  the  reaction  force  at  the  fixed  end 
and  the  time  instant  at  which  the  wave  reaches  the  fixed  end 
grows  in  magnitude  with  each  wave  reflection  at  the  fixed  end.  For 
case  #3,  the  magnitudes  of  the  peak  reaction  force  and  peak  kinetic 
energy  are  smaller  than  the  baseline  case  because  of  the  lost  yarn 
mass,  even  though  the  time  instant  of  occurrence  matches  the 
baseline  case  since  the  longitudinal  wave  speed  remains  the  same. 
Therefore  Fig.  4d  does  not  display  the  error  in  time  for  case  #3. 
The  error  in  the  predictions  of  the  reaction  force  of  case  #3  grows 
in  magnitude  with  each  reflection  at  the  fixed  end  as  can  be  seen 
from  Fig.  4c.  There  is  an  approximately  6.5%  error  in  the  reaction 
force  predictions  of  case  #4  with  respect  to  the  analytical 
predictions,  while  that  of  case  #3  is  approximately  12.2%.  This 
supports  the  conclusion  that  the  new  method  of  modeling  yarns 
with  shell  elements  is  more  accurate  compared  to  the  previous 
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-  (1 )  Baseline  with  solid  elements 

-  (2)  HEA  method  with  shell  elements 


c  £  H 
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(#3)  Force  (N)  a  (#3)  Force  (%) 

(#4)  Force  (N)  ▼  (#4)Force(%) 

(#4)  Time  (ms)  ♦  (#4)Time(%) 


#  Wave  Reflections  at  Fixed  End 


Fig.  4.  Comparison  of  results:  (a)  yarn  internal  energy,  (b)  yarn  kinetic  energy,  (c)  reaction  force  at  clamped  end,  (d)  error  in  the  predictions. 


method  with  thickness  jumps,  and  that  even  modifying  the  yarn 
material  density  in  an  attempt  to  conserve  mass  leads  to  incorrect 
results.  An  additional  advantage  is  that  the  new  method  also  better 
represents  the  yarn  cross  section.  As  we  will  see  in  Section  5,  these 
shell  elements  used  to  model  the  yarns  find  use  in  regions  away 
from  the  impact  zone  where  the  response  is  predominantly  tensile 
in  nature  and  there  are  lesser  extents  of  inter  yarn  interactions  such 
as  sliding  and  reorientation,  making  it  an  important  requirement 
that  the  tensile  response  is  accurately  captured.  We  have  just 
demonstrated  through  these  tensile  test  simulations  that  this  new 
method  of  modeling  yarns  with  shell  elements  is  capable  of  doing 
exactly  that. 

3.  Hybrid  element  analysis  applied  to  a  single  yarn 

When  choosing  between  solid  and  shell  elements,  there  is 
a  tradeoff  between  accuracy  and  computational  expense.  Therefore 
a  good  approach  would  be  to  combine  the  benefits  of  both  finite 
elements  into  the  yarn  model.  At  the  impact  zone  a  higher  level  of 
accuracy  is  needed  as  there  are  increased  projectile-curvature  and 
transverse  effects  such  as  transverse  yarn  compression,  shearing, 
and  bending  which  are  primarily  due  to  projectile-yarn  interac¬ 
tions.  At  locations  removed  from  this  impact  zone,  there  are  much 
lesser  extents  of  projectile-fabric  and  yarn-yarn  interactions.  The 
main  mode  of  yarn  deformation  is  tensile  in  nature  which  can  be 
modeled  accurately  using  a  lower  level  of  resolution.  In  this  region 
the  emphasis  is  therefore  on  computational  efficiency,  so  that  large 
yarn  dimensions  can  be  modeled.  Further,  as  seen  in  experimental 
transverse  impact  testing  of  yarns  and  fabrics,  the  predominant 
failure  is  always  underneath  the  projectile.  A  useful  approach  then 


would  be  to  use  solid  elements  to  model  the  yarns  at  the  impact 
zone  and  shell  elements  at  regions  away  from  the  impact  zone.  We 
refer  to  this  approach  as  the  hybrid  element  analysis  and  define  it 
as  ‘ the  finite  element  analysis  of  a  structure  by  combining  different 
finite  element  formulations  at  both  a  single  and  multiple  scales  of 
modeling'.  For  example,  the  yarn  is  considered  as  a  single  scale 
model  when  a  homogenized  approach  is  used  so  that  filament  level 
architecture  is  not  considered.  While  modeling  a  fabric,  multiple 
scales  of  modeling  would  imply  modeling  using  both  yarn  level 
resolution  and  a  homogenized  membrane  type  assumption  for  far 
field  regions.  In  this  study,  we  apply  the  FIEA  method  to  a  single 
scale  of  modeling  since  we  only  consider  a  yarn  level  architecture. 
Future  work  will  present  the  FIEA  method  applied  to  multiple  scales 
of  modeling.  Fig.  5  displays  the  HEA  method  applied  to  a  single 
yarn.  The  solid  elements  are  located  around  the  center  of  the  yarn 
which  is  impacted,  and  is  surrounded  on  either  side  by  shell 
elements.  A  tied  interface  is  used  between  the  two  types  of  finite 
elements.  If  an  even  number  of  solid  elements  is  used  through  the 
thickness,  it  is  also  possible  for  the  solid  and  shell  elements  to  share 
common  nodes  at  the  interface.  In  such  a  case,  the  element 
formulation  of  the  solid  elements  must  support  nodal  rotations 
which  will  ensure  compatibility  with  the  degrees  of  freedom  of  the 
shell  elements.  Since  the  material  properties  used  in  the  material 
models  for  both  the  solid  and  shell  elements  are  the  same,  and  the 
yarn  cross  sectional  area  is  the  same  by  virtue  of  our  new  modeling 
approach  for  yarns  modeled  using  shell  elements,  there  is 
a  matching  of  impedances  across  the  interface.  This  is  important  as 
it  ensures  that  there  are  no  reflections  of  the  longitudinal  strain 
wave  at  the  interface,  and  that  the  transverse  displacement  or 
bending  wave  propagates  properly  across  the  interface. 
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Fig.  5.  (a)  HEA  of  a  yarn  using  shells  and  solids  (top  view),  (b)  Close  up  at  the  interface. 


4.  Impedance  matching  in  the  HEA  approach 

Impedance  matching  is  a  crucial  requirement  especially  in  multi 
scale  models  that  contain  interfaces.  Failure  to  match  impedances 
across  interfaces  can  introduce  significant  errors  in  the  analysis  due 
to  interfacial  wave  reflections,  especially  in  wave  dominated  high 
rate  impact  phenomena  such  as  the  transverse  impact  of  woven 
fabrics.  This  was  seen  in  the  multiscale  fabric  model  of  Barauskas 
and  Abraitiene  [13,14]  where  yarn  failure  had  to  be  artificially 
‘turned  off  for  the  shell  yarn  region  adjacent  to  the  membrane 
region.  This  was  attributed  to  the  possible  overestimation  by  LS- 
DYNA  of  local  strains  in  the  finely  meshed  yarn  region  adjacent  to 
the  coarsely  meshed  membrane  region,  however  an  impedance 
mismatch  between  the  two  reasons  was  more  likely  the  primary 
cause  of  the  premature  yarn  failure,  due  a  stress  buildup  caused  by 
interfacial  reflections  of  the  longitudinal  wave,  which  then  neces¬ 
sitated  the  removal  of  the  failure  criterion.  To  demonstrate  that 
there  are  no  reflections  of  the  longitudinal  strain  wave  at  the 
shell-solid  interface  using  the  HEA  approach,  and  to  highlight  the 
detrimental  effects  produced  by  an  artificial  modeling-induced 
impedance  mismatch  at  the  interface,  four  test  cases  are  set  up 
using  the  test  case  parameters  and  yarn  material  properties  as 
before,  and  simulated  using  the  dynamic  finite  element  code  LS- 
DYNA.  In  the  first  test  case,  entitled  ‘solid’,  the  yarn  is  modeled 
using  solid  elements  and  this  serves  as  the  baseline.  In  the  second 
test  case  entitled  ‘HEA’,  the  HEA  approach  as  outlined  earlier  is 
used,  see  Fig.  5a.  The  central  portion  of  the  yarn  is  modeled  with 
solid  elements  and  measures  20  mm  while  shell  elements  are  used 
on  either  side  with  each  shell  element  region  measuring  15  mm.  In 
the  third  test  case,  entitled  ‘impedance  mismatch  #1  (IM1)’  the 
density  of  the  yarn  in  the  shell  element  region  has  been  increased 
by  1.5  times,  from  its  original  value  of  1440  kg/m3  to  2160  kg/m3.  In 
the  fourth  test  case,  entitled  ‘impedance  mismatch  #2  (IM2)’  the 
cross  sectional  area  (A)  of  the  yarn  in  the  shell  element  region  has 
been  increased  by  1.5  times,  from  its  original  value  of  0.0512  mm2 
to  0.0760  mm2.  This  dynamic  tensile  test  is  dominated  by  the 
tension  generated  in  the  yarn  due  to  the  longitudinal  elastic 
modulus.  The  inertial  effects  are  not  significant  as  evident  by  the 
yarn  kinetic  energy  which  is  negligible  compared  to  the  yarn 
internal  energy.  Further,  a  constant  velocity  is  prescribed  at  the 
pulled  end.  Thus  even  though  the  mass  of  the  yarn  changes  as  we 
are  changing  the  yarn  material  density  and  cross  sectional  area,  this 
process  allows  us  to  isolate  the  effect  the  associated  impedance 
mismatch  will  have  on  the  system  response  for  the  purpose  of 
illustration.  The  total  acoustic  impedance  is  given  by 

z  =  pcA  (13) 

In  the  third  test  case,  by  changing  the  yarn  density  of  the  shell 
element  region,  the  corresponding  longitudinal  strain  wave  speed 
will  change  according  to  Eq.  (3).  However  in  the  fourth  test  case,  the 
longitudinal  strain  wave  speed  in  both  the  solid  and  shell  element 
region  will  remain  the  same.  Fig.  6  displays  contours  of  the  longi¬ 
tudinal  tensile  stress  in  the  yarn  during  the  first  7  ps,  during  which 


time  the  wave  travels  from  the  pulled  to  the  fixed  end  and  reflects 
back  towards  the  pulled  end.  The  illustrated  labels  beside  each  yarn 
in  Fig.  6  correspond  to  ‘test  case  -  time  instant  (ps)’.  As  can  be  seen, 
the  results  of  the  HEA  case  closely  follow  the  baseline  3-d  case.  The 
longitudinal  tensile  stress  developed  behind  the  wave  front  is 
constant  as  indicated  by  the  solid  color  band  in  the  plots.  On  its  way 
to  the  fixed  end,  the  longitudinal  wave  passes  through  two 
shell-solid  interfaces,  and  as  is  observed  from  the  first  two  test 
cases,  there  are  no  wave  reflections  at  the  interfaces.  In  test  case 
IM1,  the  longitudinal  wave  lags  behind  the  other  cases  as  the  yarn 
density  has  been  increased  in  the  shell  element  region.  The  longi¬ 
tudinal  waves  for  the  cases  solid,  HEA,  and  IM2  reach  the  first 
shell-solid  interface  at  around  2  ps.  Upon  close  examination  of  the 
plots  labeled  IM1-3  and  IM1-4,  the  reflection  of  the  longitudinal 
wave  at  the  first  shell-solid  interface  is  apparent  as  the  yarn  stress 
contour  plot  no  longer  remains  a  single  solid  color  behind  the  wave 
front,  which  would  have  represented  a  constant  stress,  rather  there 
is  a  band  of  stresses  present  as  indicated  by  the  different  band  of 
colors.  Similarly,  the  plots  labeled  IM2-3  and  IM2-4  also  exhibit 
a  reflection  of  the  longitudinal  wave  at  the  interface.  Upon  reaching 
the  pulled  end,  this  reflected  longitudinal  wave  reflects  once  more 
and  heads  towards  the  fixed  end,  where  it  will  meet  the  original 
longitudinal  wave  which  is  on  its  way  back  from  the  fixed  end.  In 
addition  to  this,  there  are  also  reflections  that  will  occur  at  the 
second  shell-solid  interface.  This  leads  to  a  non-uniform  state  of 
stress  within  the  yarn  as  seen  in  the  final  plots  labeled  IM1-7  and 
IM2-7  where  a  discrete  band  of  colors  is  observed  throughout  the 
yarn  length.  This  process  continues  during  the  entire  dynamic 
tensile  test.  However  as  seen  in  the  plots  labeled  solid-7  and  HEA- 7, 
the  stress  remains  uniform  behind  the  front  of  the  longitudinal 
wave  as  indicated  by  the  single  solid  color  band.  Fig.  7  helps 
quantify  the  differences  in  the  responses  seen  in  Fig.  6.  Fig.  7a 
compares  the  longitudinal  tensile  stress  developed  at  a  location 
12.7  mm  from  the  right  pulled  end,  which  lies  within  the  shell 
element  region.  The  shape  of  the  stress  history  plot  for  both  the 
baseline  and  HEA  cases  follows  a  uniform  stepped  shape  with  the 
stress  level  remaining  constant  between  the  passing  of  the  wave 
and  the  next  arrival  of  the  reflected  longitudinal  strain  wave  from 
the  boundaries,  which  is  accompanied  by  a  sharp  incremental  rise 
in  the  stress  level.  However  as  is  evident  from  Fig.  7a,  during  the 
time  between  which  the  longitudinal  wave  has  passed  through 
a  location  and  when  it  arrives  again  at  that  location  after  a  reflec¬ 
tion  from  a  boundary,  for  both  the  IM1  and  IM2  cases,  the  stress 
level  does  not  remain  constant  and  fluctuates  with  sudden  rises 
and  dips  indicative  of  the  reflections  from  the  interfaces.  This 
causes  the  earlier  mentioned  discrete  band  of  colors  in  the  stress 
contour  plots,  seen  in  Fig.  6.  Such  an  effect  is  detrimental  as  it  can 
either  delay  or  cause  premature  yarn  failure. 

As  expected,  the  baseline  and  HEA  cases  predict  the  same  stress 
levels,  which  for  further  validation  have  been  compared  with  the 
analytical  predictions  based  on  1-d  wave  propagation  theory,  as 
outlined  earlier.  For  our  particular  test  case,  substituting  the  yarn 
material  properties  and  prescribed  velocity  at  the  pulled  end  into 
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Fig.  6.  Longitudinal  wave  propagation  through  the  yarn. 


Eq.  (10),  we  obtain  the  first  stress  increment  as  198.1  MPa.  This 
stress  value,  as  well  as  the  following  stress  increments,  closely 
agrees  with  the  baseline  and  HEA  results  with  less  than  a  percent 
difference  when  normalized  with  respect  to  the  analytical  predic¬ 
tion.  However  as  seen  from  the  IM1  and  IM2  cases  in  Fig.  7a,  the 
multiple  reflections  at  the  two  shell-solid  interfaces  cause  a  fluc¬ 
tuation  and  therefore  an  incorrect  prediction  of  the  longitudinal 
stresses.  This  error  in  predictions  grows  with  time  and  even 
exceeds  30%  at  around  18  ps  for  IM2  with  respect  to  the  analytical 
prediction.  Fig.  7b  compares  the  reaction  force  at  the  clamped  end 
of  the  yarn.  Substituting  the  yarn  material  properties  and 
prescribed  velocity  at  the  pulled  end  into  Eq.  (7),  we  obtain  the  first 
reaction  force  increment  as  20.3  N  which  closely  agrees  with  both 
the  baseline  and  HEA  cases,  again  with  less  than  a  percent  differ¬ 
ence  when  normalized  with  respect  to  the  analytical  prediction. 
However  the  IM1  and  IM2  cases  show  a  fluctuation  in  the  predic¬ 
tions,  with  the  error  once  again  increasing  with  time. 


Thus  it  becomes  important  to  match  impedances  across  any 
interfaces  present  within  the  FE  model.  As  we  will  observe  in 
Sections  6  and  8,  yarn  failure  is  implemented  using  a  stress-based 
failure  criterion  in  LS-DYNA.  Whether  the  failure  model  is  stress-  or 
strain-based,  it  becomes  all  the  more  important  that  impedances 
are  matched  across  the  interfaces  so  that  the  model  predicts  the 
correct  stress  state,  which  otherwise  could  lead  to  premature  yarn 
failure  because  of  a  stress  buildup  caused  by  multiple  reflections 
from  a  fictitious  interface. 

5.  Results  and  discussion  of  single  yarn  modeling 

Consider  the  case  of  a  0.5  g  wedge  shaped  rigid  projectile 
transversely  impacting  a  50  mm  long  uncrimped  yarn  at  its  center 
at  200  m/s.  The  projectile  moves  along  the  y-direction.  The  yarn  is 
gripped  at  both  ends.  The  material  properties  of  the  yarn  remain 
the  same  as  before.  Two  sets  of  simulations  have  been  set  up  using 
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Fig.  7.  (a)  Longitudinal  tensile  stress  at  x  =  12.7  mm  from  the  pulled  end.  (b)  Reaction  force  at  the  fixed  end. 


the  dynamic  finite  element  code  LS-DYNA.  In  the  first  set,  the  yarns 
are  modeled  using  (a)  shell  elements,  and  (b)  solid  elements.  The 
complexity  of  failure  is  not  considered  in  this  set  of  simulations  that 
focuses  on  yarn  deformations.  This  set  of  simulations  allows  us  to 
compare  the  initial  response  of  the  shell  element  case  to  the  solid 
element  case.  In  the  second  set  of  simulations,  the  yarns  are 
modeled  using  (a)  solid  elements,  and  (b)  the  HEA  approach  as 


outlined  earlier.  Failure  is  implemented  for  this  set  using  an 
element  erosion  failure  model  with  a  maximum  principal  stress 
failure  criterion  of  2759  MPa. 

Fig.  8  compares  the  results  from  the  first  set  of  simulations.  From 
Fig.  8a,  we  observe  the  velocity  agrees  in  both  cases;  however  the 
yarn  internal  energy  grows  at  a  faster  rate  from  the  moment  of 
impact  onwards  for  the  shell  element  case  and  the  yarn  kinetic 


Fig.  8.  Comparison  between  shell  elements  and  solid  elements  applied  to  transverse  yarn  impact:  (a)  projectile  velocity  and  energy  transformation  history,  (b)  projectile-yarn 
contact  force  history,  (c)  yarn  stress  history. 
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Fig.  9.  Comparison  between  solid  elements  and  HEA  applied  to  transverse  yarn  impact:  (a)  projectile  velocity  and  energy  transformation  history,  (b)  projectile-yarn  contact  force 
history,  (c)  yarn  stress  history. 


energy  begins  to  deviate  slightly  from  the  solid  element  case  after 
6  ps.  Fig.  8b  compares  the  total  projectile  to  yarn  contact  force.  This 
contact  force  is  an  indication  of  the  tension  generated  in  the  yarn. 
The  vertical  component  of  the  tension  developed  within  the  yarn  is 


Shell  Solid  HEA 

Yarn  Model 


related  to  the  decrease  in  projectile  velocity  through  the  impul¬ 
se-momentum  equation.  While  Fig.  8a  and  b  implies  similar  yarn 
tensile  forces  are  developing  in  both  cases,  it  is  not  a  sufficient 
measure  of  the  ability  of  the  shell  elements  to  be  used  for  yarn 
modeling  at  the  impact  zone.  For  this,  we  turn  to  Fig.  8c  which 
compares  the  stress  distributions  in  the  yarn  at  two  locations.  The 
longitudinal  tensile  stress  {ax)  and  transverse  shear  stress  (rxy)  have 
been  measured  near  the  bottom  surface  of  the  yarn  directly  under 
the  impact  location  of  the  projectile.  In  addition,  the  longitudinal 
tensile  stress  is  also  measured  at  a  far  field  location  that  corre¬ 
sponds  to  1  mm  from  the  right  gripped  end.  We  observe  that  the 
stress  predictions  from  the  shell  element  case  differ  from  the  solid 
element  case  at  the  impact  zone.  The  shell  element  formulation 
used  here,  as  in  most  commercial  FE  codes,  is  based  on  a  first  order 
shear  deformation  theory.  This  implies  the  transverse  shear  strains 
are  constant  through  the  thickness.  For  isotropic  materials,  a  shear 
correction  factor  of  5/6  is  used  in  an  attempt  to  correct  this  issue, 
that  leads  to  a  violation  of  the  zero  traction  condition  at  the  top  and 
bottom  surfaces  of  the  shell.  More  importantly,  the  shell  elements 
are  unable  to  account  for  any  thickness  changes  in  the  yarn  cross 
section  due  to  compression  by  the  projectile.  The  maximum 
interaction  is  seen  at  the  impact  zone  where  the  projectile 
compresses  and  shears  the  yarns  at  very  high  velocity.  In  addition 
the  curvature  of  the  projectile  induces  high  flexural  and  shear 
stresses  in  the  yarn  around  the  impact  zone.  Flere  momentum 
transfer  occurs  over  very  short  time  intervals  and  is  limited  to  the 


Fig.  10.  Computational  requirements. 
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Fig.  11.  FE  model  of  a  plain  weave  fabric  using  only  (a)  solid  elements,  (b)  shell  elements. 


region  within  the  front  of  the  longitudinal  stress  wave.  The  differ¬ 
ence  in  the  stress  predictions  at  the  impact  zone  demonstrates  the 
inability  of  the  shell  elements  to  properly  model  these  interactions, 
more  so  as  this  is  also  the  region  where  yarn  failure  would  occur 
during  an  impact  test.  However  at  far  field  regions,  the  deformation 
of  the  yarn  is  predominantly  tensile  in  nature.  We  observe  a  good 
agreement  in  the  stress  predictions  at  the  far  field  region.  This 
supports  the  reasoning  of  using  solid  elements  at  the  impact  zone 
where  higher  accuracy  is  required,  and  shell  elements  at  far  field 
regions  where  the  emphasis  is  on  computational  efficiency. 

Fig.  9  compares  the  results  from  the  second  set  of  simulations. 
The  simulation  run  time  has  been  increased  from  12  ps  to  20  ps 
with  failure  now  implemented.  From  Fig.  9a  we  observe  an  excel¬ 
lent  agreement  between  both  cases  prior  to  yarn  failure.  The  first 
instant  of  element  failure  for  the  solid  element  case  is  9.87  ps  and  is 
9.95  ps  for  the  HEA  case.  After  the  yarn  snaps,  the  yarn  internal 
energy  rapidly  drops  accompanied  by  an  increase  in  the  yarn 
kinetic  energy  as  the  yarn  flaps  freely.  The  projectile  velocity  now 
remains  constant  as  it  has  completely  penetrated  through  the  yarn. 
The  post-failure  response  displays  a  similar  trend  in  both  cases 
with  a  minor  difference  in  the  yarn  energies.  Fig.  9b  compares  the 
projectile-yarn  contact  force  and  shows  a  good  agreement 
between  both  cases.  In  Fig.  8c  we  had  observed  a  difference  in 
stress  predictions  at  the  impact  zone,  however  as  now  seen  in 
Fig.  9c,  there  is  a  good  agreement  between  the  stress  predictions  at 
both  the  impact  zone  and  far  field  regions. 

Fig.  10  compares  the  computational  requirements  between  the 
three  cases  in  terms  of  simulation  run  times  and  memory 
requirements.  It  is  important  to  consider  the  memory  requirements 
in  addition  to  the  run  times.  High  performance  clusters  such  as  the 
Beowulf  cluster  are  created  by  combining  many  similar  multi¬ 
processor  nodes  together.  However,  typically  the  memory  available 
per  node  ranges  from  2  to  8  GB.  If  the  simulation  initialization 
phase  in  LS-DYNA  requires  more  memory  than  that  available  to 
a  single  processor  or  node,  it  will  not  run  regardless  of  the  total 
memory  available  in  the  cluster.  This  problem  of  limited  memory 
available  per  node  is  not  an  issue  in  high  performance  systems  or 
workstations  that  use  a  globally  shared  type  of  memory  architec¬ 
ture  wherein  all  processors  have  access  to  all  the  memory.  Such 
systems  typically  have  a  memory  in  the  range  of  32-64  GB. 
However  the  cost  of  such  systems  usually  far  exceeds  that  of  a  high 
speed  cluster  of  similar  configuration. 

The  shell  element  case  has  the  best  computational  efficiency, 
but  for  reasons  stated  earlier,  cannot  be  used  for  regions  under¬ 
neath  the  projectile  at  the  impact  zone.  The  HEA  case  runs  almost 
40%  faster  than  the  solid  element  case.  This  unique  HEA  approach 
has  resulted  in  reproducing  the  predictions  of  the  baseline 
numerical  model  but  at  a  fraction  of  the  computational  expense. 


The  ratio  of  lengths  of  the  solid  element  portion  to  shell  element 
portion  in  this  particular  HEA  case  was  0.66,  and  by  reducing  this 
ratio,  the  computational  efficiency  can  be  further  increased. 
However  there  is  a  limit  to  the  extent  by  which  the  solid  element 
region  can  be  reduced.  For  example  in  this  particular  impact  case, 
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Fig.  12.  Comparison  between  baseline  models  for  lOOm/s  impact  velocity:  (a) 
projectile  velocity  history,  (b)  internal  energy  history.  LS-DYNA  contact  algorithm  for 
baseline  #1  (type  26),  baseline  #2  (type  a3),  baseline  #3  (type  4). 
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Fig.  13.  Setup  of  a  fabric  model  using  a  single  scale  HEA  approach:  (a)  yarns  -  solid  elements,  (b)  yarns  -  shell  elements,  (c)  projectile,  (d)  solid-shell  interface. 


from  a  series  of  simulations  we  observed  that  the  minimum  size  of 
the  solid  element  portion  must  be  at  least  two  to  three  times  the 
maximum  dimension  of  the  projectile’s  impact  face.  For  dimen¬ 
sions  less  than  this,  yarn  failure  was  observed  at  the  shell— solid 
interface  immediately  following  impact. 

6.  Single  scale  modeling  of  plain  weave  fabric  systems  and 
setup  of  the  baseline  numerical  models 

At  this  stage,  the  HEA  approach  has  been  validated  for  a  single 
yarn  model  by  comparing  it  to  both  a  baseline  yarn  FE  model  as  well 
as  analytical  predictions.  The  next  stage  is  to  extend  the  HEA 
approach  to  a  fabric  FE  model  comprised  of  these  yarns.  This  stage 
begins  with  first  setting  up  the  baseline  fabric  model  which  will  be 
used  for  comparison  to  the  HEA  fabric  model.  As  mentioned  earlier, 
a  single  scale  model  of  a  plain  weave  fabric  system  refers  to  the 
approach  where  one  level  of  resolution  is  used  throughout  the 
model.  The  fabric  may  be  entirely  modeled  by  explicitly  capturing 
the  yarn  architecture,  with  the  use  of  only  solid  elements  [19],  shell 
elements  [16],  or  1-d  elements  [21].  The  fabric  may  also  be  entirely 
modeled  by  homogenizing  the  yarns  into  a  single  and  continuous 
membrane  type  layer,  using  only  shell  elements  [3,5].  There  is  no 
mixing  of  macro  and  micro  levels  of  architecture  in  such  a  modeling 
approach.  Fig.  11a  displays  the  FE  model  of  a  single  layer  of  a  plain 
weave  fabric  modeled  using  only  solid  elements  and  maintaining 
a  yarn  level  resolution.  Fig.  lib  displays  an  equivalent  model  using 
only  shell  elements,  which  uses  our  new  method  of  non-uniform 
nodal  thicknesses  to  capture  the  yarn  geometry.  Due  to  symmetry, 
only  one  quarter  of  the  fabric  needs  to  be  modeled.  The  baseline 
fabric  numerical  model  as  seen  in  Fig.  11a  is  set  up  using  the 
preprocessor  DYNAFAB  [17]. 

Consider  the  case  of  a  0.63  g  rigid  spherical  projectile  of  diam¬ 
eter  5.55  mm  impacting  a  101.6  mm  x  50.8  mm  balanced  plain 
weave  fabric  at  the  center.  The  projectile  moves  along  the  y- 
direction.  The  fabric  is  gripped  on  all  four  sides.  The  material 
properties  of  the  yarns  used  in  the  simulation  are:  longitudinal 
elastic  modulus  ( Ex )  of  62  GPa,  transverse  elastic  moduli  (Ey)  and 
(Ez)  of  620  MPa,  and  a  maximum  principal  stress  to  failure  of 


3400  MPa.  The  yarns  are  modeled  using  single  integration  point 
solid  elements  and  the  projectile  is  modeled  using  fully  integrated 
shell  elements.  Due  to  symmetry  only  one  quarter  of  the  fabric  has 
been  modeled.  A  static  frictional  coefficient  of  0.18  is  specified  for 
the  yarn  to  yarn  contact  algorithms  and  0.18  for  the  projectile  to 
yarn  contact  algorithms.  The  thickness  of  the  warp  and  fill  yarns  is 
0.115  mm  and  the  total  fabric  thickness  at  the  cross-over  locations 
is  0.23  mm.  The  yarn  count  in  both  warp  and  fill  directions  is  34 
yarns  per  inch.  Two  impact  velocities  of  100  m/s  and  200  m/s  are 
selected. 

One  of  the  most  challenging  tasks  during  the  finite  element 
analysis  of  fabric  impact  is  to  accurately  account  for  the  highly 
complex  projectile-yarn  and  yarn-yarn  interactions  by  virtue  of 
the  complex  fabric  architecture.  This  is  accomplished  through 
contact  algorithms  in  the  FE  code  that  prevent  interpenetrations  of 
nodes  and  segments  by  applying  forces  between  penetrating  nodes 
and  segments.  Depending  on  the  contact  algorithm  used,  the 
manner  in  which  the  detection  and  removal  of  interpenetrations  is 
accomplished,  determination  of  contact  and  sliding  energies,  as 


Fig.  14.  Close  up  at  the  impact  region. 
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Fig.  15.  Comparison  of  deformation  profiles  for  the  100  m/s  impact  velocity  case  at  the  time  instants  of,  top-bottom:  (a)  15  ps,  (b)  25  ps,  (c)  40  ps,  (d)  55  ps,  (e)  60  ps,  and  (f)  70  ps. 
(Left)  Baseline  #1,  (right)  HEA. 
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Fig.  16.  Comparison  between  baseline  and  HEA  models  for  impact  velocity  of  100  m/s:  (a)  projectile  velocity,  (b)  fabric  internal  energy,  (c)  sliding  energy,  (d)  contact  force.  LS-DYNA 
contact  algorithm  for  baseline  #1  (type  26),  baseline  #2  (type  a3). 


well  as  the  transfer  of  boundaries  between  interacting  segments 
vary.  Consequently  with  the  use  of  different  contact  algorithms,  the 
system  response  and  computational  run  times  for  the  same  phys¬ 
ical  problem  may  differ  either  negligibly  or  significantly.  It  is 
therefore  important  to  investigate  this  effect  and  accordingly  three 
baseline  simulations  are  set  up  using  LS-DYNA  contact  types  26,  a3, 
and  4  respectively.  Further  details  of  the  implementation  and 
working  of  these  contact  algorithms  can  be  obtained  from  Ref.  [22]. 
The  first  baseline  uses  the  most  robust  contact  algorithm  of  the 
three  and  consequently  is  the  most  computationally  intensive, 
while  the  second  baseline  uses  the  least  computationally  intensive 
algorithm  for  this  type  of  problem.  The  first  and  third  baselines  use 
a  single  surface  type  contact  algorithm  which  ensures  that  no  two 
adjacent  yarns  (viz.  warp-warp  or  fill-fill)  can  penetrate  each 
other.  The  second  baseline  uses  a  master  and  slave  surface  type 
algorithm  which  means  that  the  warp  yarns  will  not  penetrate  the 
fill  yarns  and  vice-versa,  but  there  may  be  a  possibility  of  two  warp 
or  two  fill  yarns  penetrating  each  other  should  they  ever  come  into 
contact.  For  these  single  layer  impact  simulations,  the  first  baseline 
serves  to  be  the  most  robust.  Flowever  as  the  size  of  this  fabric 
model  becomes  larger,  the  first  baseline  becomes  computationally 
very  intensive  compared  to  the  other  baselines  in  terms  of  both 
memory  requirements  and  processing  power.  Also  the  third  base¬ 
line  remains  intensive  in  terms  of  memory  even  though  the  pro¬ 
cessing  power  requirements  compare  to  that  of  the  second 
baseline.  In  such  cases,  the  less  intensive  contact  algorithm  of  the 
second  baseline  becomes  preferable.  These  differences  in 


computing  requirements  between  the  three  baseline  cases  have 
been  quantified  and  presented  in  Section  8.  It  is  important  to  note 
here  that  with  multiple  fabric  layers,  an  eroding  type  of  contact  (LS- 
DYNA  contact  types  14  or  15)  becomes  necessary  to  prevent 
penetrations  between  the  projectile  and  the  fabric  as  well  as 
between  the  various  fabric  layers,  while  the  projectile  is  pene¬ 
trating  through.  This  type  of  contact  algorithm  is  even  more 
intensive  than  the  first  baseline  used  here.  Flowever  by  using 
different  contact  algorithms  in  different  regions,  much  like  our  FIEA 
model  uses  different  finite  elements  and  modeling  resolutions  in 
different  regions  of  the  fabric,  the  overall  computational  require¬ 
ments  can  be  reduced.  One  example  would  be  to  include  a  single 
surface  eroding  type  contact  algorithm  for  the  region  around  the 
impact  zone  only,  with  a  less  intensive  type  of  contact  algorithm  at 
far  field  regions  from  the  impact  zone  where  element  erosion  and 
penetration  between  two  adjacent  yarns  is  unlikely  to  occur.  There 
is  no  clear  choice  as  to  which  contact  algorithm  can  be  considered 
the  best  to  use  for  all  cases,  since  each  contact  algorithm  has  its 
own  advantages  as  stated  earlier.  The  choice  will  depend  on  the 
parameters  and  physics  of  the  impact  problem  simulated  and 
finally  this  is  left  to  the  experience  of  the  FE  analyst. 

Fig.  12  compares  the  projectile  velocity  and  fabric  internal 
energy  histories  between  the  three  baseline  models  for  the  100  m/s 
impact  velocity  case.  The  projectile  is  able  to  completely  penetrate 
the  fabric  target  with  both  impact  velocity  cases.  The  first  instant  of 
element  erosion  is  slightly  different  in  each  baseline  which  triggers 
subsequent  element  failure  at  different  times.  Thus  the  projectile 
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Fig.  17.  Comparison  between  baseline  and  HEA  models  for  impact  velocity  of  200  m/s: 
(a)  projectile  velocity,  (b)  fabric  internal  energy.  LS-DYNA  contact  algorithm  for 
baseline  #1  (type  26),  baseline  #2  (type  a3). 


velocity  and  fabric  internal  energy  histories  until  the  first  instance 
of  failure  are  exactly  the  same  in  all  three  baselines.  However,  even 
after  complete  penetration,  the  residual  velocity  of  the  projectile 
closely  agrees  between  all  baselines.  The  deformation  of  the  fabric 
is  pyramidal  in  nature,  as  seen  in  subsequent  contour  plots  of  fabric 
deformation  profiles,  and  this  agrees  with  experimental  observa¬ 
tions.  The  shape  of  the  base  of  this  pyramid  is  obtained  by  joining 
the  front  of  the  transverse  displacement  wave  in  each  yarn  by  a  line. 
When  the  outward  propagating  longitudinal  strain  wave  reaches 
the  four  clamped  boundaries,  it  reflects  back  towards  the  projectile. 
With  each  interaction  of  the  longitudinal  wave  and  the  transverse 
displacement  wave,  the  velocity  of  the  transverse  displacement 
wave  is  increased.  The  deformation  pyramid  does  not  reach  the 
clamped  boundaries  before  the  projectile  penetrates  through  the 
fabric  in  both  impact  velocity  cases.  While  the  projectile  is  pene¬ 
trating  through  a  hole  in  the  fabric  which  is  smaller  than  the 
projectile’s  diameter,  the  frictional  energy  dissipated  by  the  pro¬ 
jectile-yarn  interaction  rises  to  its  maximum  level.  While  the  fabric 
is  springing  back  after  complete  projectile  penetration,  the  fric¬ 
tional  energy  dissipated  by  yarn-yarn  sliding  interactions  rises  to 
its  maximum  level. 

This  fabric  model  where  the  explicit  yarn  architecture  is 
captured  with  solid  elements  is  used  as  a  baseline  numerical  model 


with  which  to  compare  the  HEA  models  from  subsequent  sections. 
The  HEA  models  will  be  compared  to  both  the  first  and  second 
baselines  which  respectively  correspond  to  the  most  robust  and 
computationally  intensive  (LS-DYNA  type  26),  and  least  computa¬ 
tionally  intensive  contact  algorithm  (LS-DYNA  type  a3).  The  focus 
will  be  to  reproduce  the  baseline  fabric  system  response  but  at 
a  lower  computational  expense. 

7.  Hybrid  element  analysis  applied  to  a  fabric 

In  multi  scale  modeling  with  the  HEA  method,  different  reso¬ 
lution  levels  modeled  with  different  finite  elements  are  incorpo¬ 
rated  into  one  single  model.  In  this  paper,  the  hybrid  element 
analysis  has  first  been  applied  to  create  single  scale  models,  where 
the  fabric  is  modeled  using  a  yarn  level  architecture  with  both  solid 
and  shell  elements.  This  combines  the  accuracy  of  yarn  level 
resolution  with  solid  elements,  with  the  computational  efficiency 
of  shell  elements.  The  region  modeled  with  a  yarn  level  architec¬ 
ture  is  referred  to  as  the  local  region.  Future  work  will  deal  with 
multi  scale  HEA  models.  In  the  future  model,  a  homogenized 
membrane  type  region  will  be  added  to  the  model  at  far  field 
locations  and  is  referred  to  as  the  global  region.  Fig.  13  illustrates  the 
various  regions  of  a  typical  single  scale  HEA  model.  An  explicit  yarn 
architecture  is  used  everywhere.  For  reasons  cited  earlier,  solid 
elements  have  been  used  around  the  impact  zone  while  shell 
elements  have  been  used  everywhere  else.  The  HEA  approach 
always  ensures  that  the  propagation  of  the  longitudinal  strain  wave 
and  transverse  displacement  wave  across  the  solid-shell  yarn 
interface  is  unaffected  due  to  the  impedance  matching  across  the 
solid-shell  interface.  Fig.  14  is  a  close  up  at  the  impact  zone,  and 
also  displays  the  mesh  scheme.  Due  to  symmetry  only  one  quarter 
of  the  fabric  has  been  modeled. 

8.  Results  and  discussion 

A  single  scale  HEA  model  is  set  up  to  compare  results  against  the 
baseline  numerical  models  described  in  Section  6.  Consider  the 
same  fabric  impact  cases  as  the  baselines,  where  the  yarns  have  the 
same  material  properties.  Yarn  failure  has  been  incorporated  using 
an  element  erosion  model  with  a  maximum  principal  stress  failure 


Y///A  Memory  Requirements 


i  1  i  •  i  i 

Baseline  #1  Baseline  #2  Baseline  #3  Single  Scale  HEA 

Fabric  Model 


Fig.  18.  Simulation  memory  requirements  and  run  times. 
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criterion  of  3300  MPa.  For  the  HEA  model,  the  total  dimensions  of 
the  central  square  patch  with  local  yarn  architecture  modeled  using 
solid  elements  are  9.5  mm  x  9.5  mm  while  the  remaining  fabric 
area  consists  of  local  yarn  architecture  modeled  using  shell 
elements.  The  ratio  of  the  areas  of  the  local  solid  element  region  to 
the  local  shell  element  region  is  0.018,  and  that  of  the  shell  element 
region  to  total  fabric  area  is  0.982.  LS-DYNA  contact  type  26  is  used 
for  contact  between  the  solid  element  yarns  while  contact  type  a3 
is  used  for  the  shell  element  yarns.  The  interface  between  the  solid 
and  shell  element  yarns  is  created  using  the  LS-DYNA  keyword 
*CONSTRAINED_SHELL_TO_SOLID.  The  preprocessor  DYNA-HEA  [17] 
is  used  to  automatically  create  the  FE  mesh  and  all  interface  defi¬ 
nitions  of  the  FIEA  model. 

Fig.  15  compares  the  contours  of  vertical  displacement  between 
the  first  baseline  and  the  HEA  cases,  as  seen  from  the  top  view  for 
the  100  m/s  impact  velocity  case.  These  contours  are  indicative  of 
the  propagation  of  the  transverse  displacement  wave.  The  pyra¬ 
midal  deformation  profile  agrees  with  experimental  observations. 
There  is  a  good  agreement  between  the  deformation  profiles.  The 
transverse  wave  propagates  slightly  faster  towards  the  two  longer 
sides  of  the  fabric  because  of  the  larger  number  of  reflections  of  the 
longitudinal  wave  at  the  longer  boundaries  compared  to  the 
shorter  boundaries.  This  implies  that  the  number  of  interactions 
between  the  longitudinal  wave  and  transverse  displacement  wave 
is  greater  along  the  direction  of  the  shorter  sides.  The  speed  of  the 
transverse  displacement  wave  is  increased  with  each  interaction. 
The  transverse  displacement  wave  crosses  the  local-global  inter¬ 
face  by  the  time  the  projectile  has  penetrated  through  the  fabric. 

Fig.  16  compares  the  time  history  results  of  the  baseline  and  HEA 
cases  for  the  100  m/s  impact  velocity  case.  There  is  a  very  good 
agreement  between  the  projectile  velocity  and  fabric  internal 
energy  responses,  especially  prior  to  failure  as  seen  in  Fig.  16a  and 
b.  Fig.  16c  compares  the  total  sliding  energy  which  is  the  energy 
dissipated  when  two  surfaces  slide  past  one  another.  This  includes 
both  projectile-yarn  and  yarn-yarn  interactions.  Since  the  fabric  is 
gripped  on  all  four  sides,  the  sliding  energy  remains  small  until  the 
projectile  begins  to  penetrate  through  the  fabric.  After  complete 
penetration,  the  fabric  begins  to  elastically  recover  or  spring  back 
towards  its  initial  shape  which  causes  the  sliding  energy  to  rapidly 
increase.  Fig.  16d  compares  the  projectile  to  fabric  contact  force  for 


all  cases.  This  provides  a  measure  of  the  resistance  force  that 
decelerates  the  projectile  as  well  as  the  tension  developed  within 
the  deformation  pyramid.  There  is  a  very  good  agreement  between 
the  predictions. 

Fig.  17  compares  time  history  results  of  the  baseline  and  HEA 
cases  for  the  200  m/s  impact  velocity  case.  The  first  as  well  as 
subsequent  instances  of  element  failure  vary  between  the  baselines 
and  the  difference  in  the  projectile  velocity  history  after  the  onset 
of  yarn  failure  becomes  apparent  starting  at  around  15  ps.  This 
small  difference  is  presumably  due  to  the  different  contact  algo¬ 
rithms  used  since  all  other  parameters  have  been  kept  the  same. 
There  is  a  good  agreement  in  the  projectile  velocity  and  fabric 
internal  energy  histories  especially  between  the  second  baseline 
and  the  HEA  cases. 

Fig.  18  compares  the  computational  requirements  for  both 
impact  velocity  cases.  This  consists  of  both  the  total  run  time  as  well 
as  the  memory  requirements  of  all  models.  The  simulations  were 
run  using  LS-DYNA  SMP  version  971  on  a  64-bit  Dell  Precision  690 
workstation  with  four  Intel  Xeon  3.00  GHz  processors  and  16  GB  of 
available  RAM.  The  single  scale  HEA  model  with  fully  local  yarn 
architecture  ran  faster  than  all  three  baseline  cases,  and  was 
approximately  2  times  faster  than  the  first  baseline.  These  savings 
are  not  just  reflected  in  the  run  times  but  also  in  the  memory 
requirements  of  each  model  as  seen  in  Fig.  18,  where  the  HEA 
model  required  the  least  amount  of  memory.  This  is  also  an 
important  consideration  for  reasons  cited  earlier. 

This  savings  in  computational  expense  while  still  reproducing 
the  system  response  of  the  baseline  simulations  demonstrates  the 
usefulness  of  the  HEA  approach.  As  the  dimensions  of  the  model 
and  the  simulation  termination  time  increase,  these  savings 
become  even  larger.  Fig.  19  displays  the  percentage  utilization  of 
the  central  processing  unit  (CPU)  of  the  computer  with  respect  to 
the  run  times  reported  in  Fig.  18.  In  the  first  baseline  case,  we 
observe  that  the  contact  algorithm  utilized  more  CPU  resources 
than  the  processing  of  the  finite  elements.  This  was  expected  since 
this  case  corresponds  to  the  most  robust  contact  algorithm.  We  see 
this  trend  reversed  in  the  other  two  baseline  cases  that  used 
computationally  less  intensive  contact  algorithms  while  keeping 
the  total  number  of  finite  elements  and  their  formulation  the  same. 
Of  all  the  cases,  we  observe  that  percentage  CPU  utilization  of  the 
element  processing  outweighed  that  of  the  contact  algorithms  the 
greatest  in  the  single  scale  HEA  case.  Further  even  though  this  case 
used  a  fully  local  architecture,  the  relative  CPU  utilization  by  the 
contact  algorithm  was  far  less  than  that  of  the  baselines,  because 
the  contact  pairs  defined  were  predominantly  for  shell  yarn  to  shell 
yarn  unlike  the  solid  yarn  to  solid  yarn  pairs  in  the  baseline  cases. 

9.  Conclusions 

The  finite  element  modeling  of  the  impact  of  flexible  plain 
weave  fabrics  was  discussed.  A  baseline  fabric  model  comprising 
a  yarn  level  architecture  modeled  with  solid  elements  was  pre¬ 
sented.  The  use  of  different  contact  algorithms  between  the  yarns 
in  the  finite  element  model  was  shown  to  have  an  effect  on  the 
fabric  system  response  as  well  as  the  computational  requirements. 
A  novel  approach,  the  hybrid  element  analysis,  was  introduced  that 
incorporated  using  different  finite  element  formulations  currently 
at  a  single  scale  of  modeling.  Solid  elements  were  used  to  model 
yarns  at  the  impact  zone  while  shell  elements  were  used  else¬ 
where.  This  choice  of  the  appropriate  finite  element  formulations 
to  use  in  the  various  regions  of  the  fabric  model  was  based  on  the 
demonstrated  inability  of  the  shell  elements  to  capture  the  correct 
stress  response  and  yarn  thickness  changes  in  the  impact  zone 
directly  underneath  the  projectile.  The  detrimental  effects  caused 
by  a  mismatch  of  impedance  at  the  solid-shell  interface  were 
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highlighted.  The  HEA  approach  was  able  to  accurately  reproduce 
the  baseline  simulation  results  but  at  a  fraction  of  the  computa¬ 
tional  expense,  proving  it  to  be  an  invaluable  tool  in  the  modeling 
and  simulation  of  the  impact  of  fabric  systems.  This  research  paper 
dealt  with  single  scale  modeling  issues.  Future  work  will  deal  with 
multi  scale  modeling  issues  using  the  HEA  method,  wherein  an 
even  greater  savings  in  computational  expense  is  envisioned  with 
the  inclusion  of  the  homogenized  or  global  region. 
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